##########################################################################################

library(data.table)
library(optparse)
library(ArchR)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/subcell/cluster2/cluster2.combineRNA.motif_peak2gene.Rdata"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot/sperm_enhancer/"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

a <- load(comine_data_file)


corrCutoff <- 0.5       # Default in plotPeak2GeneHeatmap is 0.45
varCutoffATAC <- 0.25   # Default in plotPeak2GeneHeatmap is 0.25
varCutoffRNA <- 0.25    # Default in plotPeak2GeneHeatmap is 0.25

###########################################################################################
nclust <- 25
atac_proj <- addPeak2GeneLinks(ArchRProj = atac_proj , useMatrix = "GeneExpressionMatrix")

p2gMat <- plotPeak2GeneHeatmap(
  atac_proj, 
  corCutOff = corrCutoff, 
  groupBy="cell_type",
  nPlot = 1000000, returnMatrices=TRUE, 
  k=nclust, seed=1)

## 输出peak-gene相关参数
kclust_df_out <- cbind( kclust=p2gMat$ATAC$kmeansId , p2gMat$Peak2GeneLinks )
tmp_peak_dat <- data.frame(atac_proj@peakSet)
tmp_peak_dat$peak <- paste0( tmp_peak_dat$seqnames , ":" , tmp_peak_dat$start , "-" , tmp_peak_dat$end )
kclust_df_out <- merge( kclust_df_out , tmp_peak_dat , by = "peak" )

out_file <- paste0(out_path, sprintf("/peakToGeneHeatmap_LabelClust_k%s.tsv", nclust))
write.table(kclust_df_out , out_file , row.names = F , sep = "\t")